My Attempt at Your Model
Brother Saunders
library(car)
library(mosaic)
library(DT)
library(tidyverse)Read in the Data
files <- dir()
rdat <- read.csv(grep(".csv", files, value=TRUE), header=TRUE)
as_tibble(rdat)# A tibble: 2,600 × 11
y x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
<dbl> <int> <int> <int> <int> <int> <int> <int> <dbl> <int> <int>
1 -4.38 0 0 0 0 0 0 0 -0.0795 1 -1
2 7.66 1 1 0 0 0 1 0 -2.54 1 -1
3 -6.82 1 1 0 0 0 1 0 -2.36 -1 -1
4 -1.23 1 1 0 0 0 1 0 1.45 -1 -1
5 4.01 0 0 0 0 0 0 0 0.690 1 1
6 0.730 0 0 0 0 0 0 0 0.175 1 -1
7 -6.80 1 1 0 0 0 1 0 1.72 1 1
8 4.99 0 0 0 0 0 0 0 -0.419 -1 1
9 5.70 1 1 0 0 0 1 0 -0.199 1 -1
10 -1.98 0 0 0 0 0 0 0 0.302 -1 -1
# ℹ 2,590 more rowsCreate first Pairs Plot
pairs(rdat, pch=16, cex=1, panel=panel.smooth, col=rgb(.8,.8,.8))
Start with everything except x1 as x1= x2
rdat <- rdat[,-2]
lm1 <- lm(y ~ ., data=rdat)
summary(lm1)##
## Call:
## lm(formula = y ~ ., data = rdat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.189 -4.180 0.037 4.136 32.647
##
## Coefficients: (4 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.23119 0.24247 -0.953 0.34044
## x2 0.17530 0.35238 0.497 0.61890
## x3 1.33716 0.46347 2.885 0.00395 **
## x4 NA NA NA NA
## x5 NA NA NA NA
## x6 NA NA NA NA
## x7 NA NA NA NA
## x8 0.31343 0.13045 2.403 0.01634 *
## x9 0.07266 0.16075 0.452 0.65131
## x10 0.14065 0.16086 0.874 0.38201
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.192 on 2594 degrees of freedom
## Multiple R-squared: 0.00613, Adjusted R-squared: 0.004214
## F-statistic: 3.2 on 5 and 2594 DF, p-value: 0.006962Sample size is large, try for more.
lm2 <- lm(y ~ .^5, data=rdat)
summary(lm2)$coefficients %>%
as.data.frame() %>%
arrange(`Pr(>|t|)`) %>%
datatable()Pull out significant terms
lm3 <- lm(y ~ x8 + x3: x8 + x2: x8: x9: x10 + x8: x9: x10 + x3: x8: x9: x10, data=rdat)
summary(lm3)##
## Call:
## lm(formula = y ~ x8 + x3:x8 + x2:x8:x9:x10 + x8:x9:x10 + x3:x8:x9:x10,
## data = rdat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.8110 -4.1093 -0.0493 4.2569 25.1425
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.08016 0.12479 -0.642 0.52071
## x8 0.31071 0.11487 2.705 0.00688 **
## x8:x3 -0.49795 0.24344 -2.045 0.04091 *
## x8:x9:x10 -6.89563 0.23465 -29.387 < 2e-16 ***
## x8:x2:x9:x10 7.22481 0.26916 26.842 < 2e-16 ***
## x8:x3:x9:x10 13.20240 0.31808 41.507 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.355 on 2594 degrees of freedom
## Multiple R-squared: 0.402, Adjusted R-squared: 0.4009
## F-statistic: 348.8 on 5 and 2594 DF, p-value: < 2.2e-16palette(c("skyblue","orange","green","purple","steelblue","red","green3","black"))
pairs(cbind(R=lm3$res, Fit=lm3$fit, rdat), pch=16, cex=1, panel=panel.smooth, col=interaction(rdat$x2, rdat$x3, rdat$x9, rdat$x10, drop=TRUE))
Zoom in
ggplot(rdat, aes( x= x8, y=y, color=interaction( x2, x3, x9, x10))) +
geom_point() +
facet_wrap(~interaction( x2, x3, x9, x10))
ggplot(rdat, aes( x= x8, y=y, color=interaction( x2, x3, x9, x10, x4))) +
geom_point() +
facet_wrap(~interaction( x2, x3, x9, x10))
ggplot(rdat, aes( x= x8, y=y, color=interaction( x2, x3, x9, x10, x5))) +
geom_point() +
facet_wrap(~interaction( x2, x3, x9, x10))
ggplot(rdat, aes( x= x8, y=y, color=interaction( x2, x3, x9, x10, x6))) +
geom_point() +
facet_wrap(~interaction( x2, x3, x9, x10))
ggplot(rdat, aes( x= x8, y=y, color=interaction( x2, x3, x9, x10, x7))) +
geom_point() +
facet_wrap(~interaction( x2, x3, x9, x10))
ggplot(rdat, aes( x= x8, y=y, color=interaction( x2, x3, x9, x10))) +
geom_point() +
facet_wrap(~interaction( x2, x3, x9, x10))
Try my own switches
rdat2 <- rdat %>%
mutate(mycolor = case_when(
x2==0 & x3==0 & x9==-1 & x10==-1 ~ "Cubic Down",
x2==0 & x3==0 & x9==1 & x10==1 ~ "Cubic Down",
x2==0 & x3==0 & x9==1 & x10==-1 ~ "Cubic Up",
x2==0 & x3==0 & x9==-1 & x10==1 ~ "Cubic Up",
x2==1 & x3==0 & x9==-1 & x10==-1 ~ "Quintic Up",
x2==1 & x3==0 & x9==1 & x10==1 ~ "Quintic Up",
x2==1 & x3==0 & x9==1 & x10==-1 ~ "Quintic Down",
x2==1 & x3==0 & x9==-1 & x10==1 ~ "Quintic Down",
x2==0 & x3==1 & x9==-1 & x10==-1 & x8>0 ~ "Quadratic Up Right",
x2==0 & x3==1 & x9==1 & x10==1 & x8>0~ "Quadratic Up Right",
x2==0 & x3==1 & x9==1 & x10==-1 & x8<0 ~ "Quadratic Up Left",
x2==0 & x3==1 & x9==-1 & x10==1 & x8<0~ "Quadratic Up Left",
x2==0 & x3==1 & x9==-1 & x10==1 & x8>0~ "Quadratic Down Right",
x2==0 & x3==1 & x9==1 & x10==-1 & x8>0~ "Quadratic Down Right",
x2==0 & x3==1 & x9==-1 & x10==-1 & x8<0~ "Quadratic Down Left",
x2==0 & x3==1 & x9==1 & x10==1 & x8<0~ "Quadratic Down Left",
))
ggplot(rdat2, aes( x= x8, y=y, color=mycolor)) +
geom_point() +
facet_wrap(~interaction( x2, x3, x9, x10))
ggplot(rdat2, aes( x= x8, y=y, color=mycolor)) +
geom_point() +
facet_wrap(~mycolor)
lm4 <- lm(y ~ mycolor + x8:mycolor + I( x8^2):mycolor + I( x8^3):mycolor + I( x8^4):mycolor + I( x8^5):mycolor, data=rdat2)
summary(lm4)##
## Call:
## lm(formula = y ~ mycolor + x8:mycolor + I(x8^2):mycolor + I(x8^3):mycolor +
## I(x8^4):mycolor + I(x8^5):mycolor, data = rdat2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.2357 -2.5761 -0.0648 2.6062 13.4635
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.14466 0.32037 -0.452 0.65164
## mycolorCubic Up 0.05787 0.44283 0.131 0.89603
## mycolorQuadratic Down Left -36.98867 24.74862 -1.495 0.13515
## mycolorQuadratic Down Right -70.42413 16.46831 -4.276 1.97e-05 ***
## mycolorQuadratic Up Left -15.26175 24.37863 -0.626 0.53135
## mycolorQuadratic Up Right 53.60584 17.47839 3.067 0.00219 **
## mycolorQuintic Down 0.04223 0.45693 0.092 0.92636
## mycolorQuintic Up -0.51824 0.45716 -1.134 0.25706
## mycolorCubic Down:x8 11.58827 0.93520 12.391 < 2e-16 ***
## mycolorCubic Up:x8 -9.30418 0.92190 -10.092 < 2e-16 ***
## mycolorQuadratic Down Left:x8 -28.64204 130.31594 -0.220 0.82605
## mycolorQuadratic Down Right:x8 228.28264 93.11730 2.452 0.01429 *
## mycolorQuadratic Up Left:x8 -190.93491 123.56317 -1.545 0.12241
## mycolorQuadratic Up Right:x8 -126.48696 94.55749 -1.338 0.18112
## mycolorQuintic Down:x8 -2.68016 0.50609 -5.296 1.29e-07 ***
## mycolorQuintic Up:x8 3.12033 0.50134 6.224 5.65e-10 ***
## mycolorCubic Down:I(x8^2) -0.24955 1.06256 -0.235 0.81434
## mycolorCubic Up:I(x8^2) 0.13598 1.00547 0.135 0.89243
## mycolorQuadratic Down Left:I(x8^2) 79.11401 248.27532 0.319 0.75001
## mycolorQuadratic Down Right:I(x8^2) -317.74366 187.34056 -1.696 0.08999 .
## mycolorQuadratic Up Left:I(x8^2) -398.50591 228.15814 -1.747 0.08082 .
## mycolorQuadratic Up Right:I(x8^2) 116.86690 183.32486 0.637 0.52387
## mycolorQuintic Down:I(x8^2) 0.08022 0.31566 0.254 0.79940
## mycolorQuintic Up:I(x8^2) 0.42432 0.29088 1.459 0.14477
## mycolorCubic Down:I(x8^3) -18.30246 1.97640 -9.261 < 2e-16 ***
## mycolorCubic Up:I(x8^3) 13.94968 1.90453 7.324 3.20e-13 ***
## mycolorQuadratic Down Left:I(x8^3) 128.44133 217.91564 0.589 0.55564
## mycolorQuadratic Down Right:I(x8^3) 229.56497 171.47661 1.339 0.18077
## mycolorQuadratic Up Left:I(x8^3) -334.51958 194.65175 -1.719 0.08582 .
## mycolorQuadratic Up Right:I(x8^3) -48.29833 161.95525 -0.298 0.76556
## mycolorQuintic Down:I(x8^3) 2.88979 0.30231 9.559 < 2e-16 ***
## mycolorQuintic Up:I(x8^3) -2.88273 0.28754 -10.025 < 2e-16 ***
## mycolorCubic Down:I(x8^4) 0.32585 0.64566 0.505 0.61383
## mycolorCubic Up:I(x8^4) -0.08951 0.59065 -0.152 0.87956
## mycolorQuadratic Down Left:I(x8^4) 65.96118 89.51113 0.737 0.46125
## mycolorQuadratic Down Right:I(x8^4) -83.82155 72.65647 -1.154 0.24874
## mycolorQuadratic Up Left:I(x8^4) -126.24773 77.80473 -1.623 0.10479
## mycolorQuadratic Up Right:I(x8^4) 7.53986 66.27789 0.114 0.90944
## mycolorQuintic Down:I(x8^4) 0.01270 0.05255 0.242 0.80901
## mycolorQuintic Up:I(x8^4) -0.05789 0.04616 -1.254 0.20988
## mycolorCubic Down:I(x8^5) 1.57592 0.93894 1.678 0.09339 .
## mycolorCubic Up:I(x8^5) 0.34243 0.87828 0.390 0.69665
## mycolorQuadratic Down Left:I(x8^5) 11.79943 13.92928 0.847 0.39702
## mycolorQuadratic Down Right:I(x8^5) 11.71172 11.55545 1.014 0.31091
## mycolorQuadratic Up Left:I(x8^5) -18.31457 11.78456 -1.554 0.12028
## mycolorQuadratic Up Right:I(x8^5) 0.35124 10.19107 0.034 0.97251
## mycolorQuintic Down:I(x8^5) -0.46370 0.03995 -11.607 < 2e-16 ***
## mycolorQuintic Up:I(x8^5) 0.45323 0.03624 12.507 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.922 on 2552 degrees of freedom
## Multiple R-squared: 0.7759, Adjusted R-squared: 0.7718
## F-statistic: 188 on 47 and 2552 DF, p-value: < 2.2e-16ggplot(rdat2, aes( x= x8, y=y, color=mycolor)) +
geom_point(pch=1) +
geom_point(aes(y=lm4$fit), cex=0.5) +
facet_wrap(~mycolor)
Reduce
Played with this for quite a while. There must be a beautiful combination of x2, x3, x8, x9 that made these all happen with your original switches. But it was easier just to make my own switches.
Well done.
rdat2 <- rdat2 %>%
mutate( x11 = ifelse(mycolor=="Quadratic Down Right", 1, 0),
x12 = ifelse(mycolor=="Quadratic Up Right", 1, 0),
x13 = ifelse(mycolor=="Cubic Down", 1, 0),
x14 = ifelse(mycolor=="Cubic Up", 1, 0),
x15 = ifelse(mycolor=="Quintic Down", 1, 0),
x16 = ifelse(mycolor=="Quintic Up", 1, 0),
x17 = ifelse(mycolor=="Quadratic Down Left", 1, 0),
x18 = ifelse(mycolor=="Quadratic Up Left", 1, 0))
lm5 <- lm(y ~ x11 + x11: x8 + x11:I( x8^2) + #Quadratic Down Right
x12 + x12: x8 + x12:I( x8^2) + #Quadratic Up Right
x17 + x17: x8 + x17:I( x8^2) + #Quadratic Down Left
x18 + x18: x8 + x18:I( x8^2) + #Quadratic Up Left
x13: x8 + x13:I( x8^3) + #Cubic Down
x14: x8 + x14:I( x8^3) + #Cubic Up
x15: x8 + x15:I( x8^3) + x15:I( x8^5) + #Quintic Down
x16: x8 + x16:I( x8^3) + x16:I( x8^5) #Quintic Up
, data=rdat2)
summary(lm5)##
## Call:
## lm(formula = y ~ x11 + x11:x8 + x11:I(x8^2) + x12 + x12:x8 +
## x12:I(x8^2) + x17 + x17:x8 + x17:I(x8^2) + x18 + x18:x8 +
## x18:I(x8^2) + x13:x8 + x13:I(x8^3) + x14:x8 + x14:I(x8^3) +
## x15:x8 + x15:I(x8^3) + x15:I(x8^5) + x16:x8 + x16:I(x8^3) +
## x16:I(x8^5), data = rdat2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.0607 -2.6280 -0.0454 2.6323 13.0411
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.05800 0.08474 -0.684 0.494
## x11 -36.12331 1.72245 -20.972 < 2e-16 ***
## x12 35.41434 1.89730 18.666 < 2e-16 ***
## x17 -34.54107 2.67490 -12.913 < 2e-16 ***
## x18 37.99485 2.35524 16.132 < 2e-16 ***
## x11:x8 56.23008 3.19009 17.626 < 2e-16 ***
## x11:I(x8^2) -22.77164 1.26613 -17.985 < 2e-16 ***
## x8:x12 -52.77655 3.24250 -16.277 < 2e-16 ***
## I(x8^2):x12 21.06233 1.21148 17.386 < 2e-16 ***
## x8:x17 -49.97923 4.33837 -11.520 < 2e-16 ***
## I(x8^2):x17 -19.87236 1.61987 -12.268 < 2e-16 ***
## x8:x18 57.71141 3.81713 15.119 < 2e-16 ***
## I(x8^2):x18 23.23116 1.38806 16.736 < 2e-16 ***
## x8:x13 10.34734 0.52746 19.617 < 2e-16 ***
## x13:I(x8^3) -15.12264 0.43847 -34.490 < 2e-16 ***
## x8:x14 -9.60837 0.51705 -18.583 < 2e-16 ***
## I(x8^3):x14 14.68082 0.41038 35.774 < 2e-16 ***
## x8:x15 -2.64218 0.50620 -5.220 1.94e-07 ***
## I(x8^3):x15 2.85507 0.30207 9.452 < 2e-16 ***
## x15:I(x8^5) -0.45953 0.03993 -11.509 < 2e-16 ***
## x8:x16 3.19876 0.50060 6.390 1.96e-10 ***
## I(x8^3):x16 -2.92284 0.28742 -10.169 < 2e-16 ***
## I(x8^5):x16 0.45776 0.03626 12.625 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.934 on 2577 degrees of freedom
## Multiple R-squared: 0.7723, Adjusted R-squared: 0.7704
## F-statistic: 397.3 on 22 and 2577 DF, p-value: < 2.2e-16ggplot(rdat2, aes( x= x8, y=y, color=mycolor)) +
geom_point(pch=1) +
geom_point(aes(y=lm5$fit), cex=0.5) +
facet_wrap(~mycolor)
Final Guess
lm
rdat2 <- rdat %>%
mutate(mycolor = case_when(
x2==0 & x3==0 & x9==-1 & x10==-1 ~ "Cubic Down",
x2==0 & x3==0 & x9==1 & x10==1 ~ "Cubic Down",
x2==0 & x3==0 & x9==1 & x10==-1 ~ "Cubic Up",
x2==0 & x3==0 & x9==-1 & x10==1 ~ "Cubic Up",
x2==1 & x3==0 & x9==-1 & x10==-1 ~ "Quintic Up",
x2==1 & x3==0 & x9==1 & x10==1 ~ "Quintic Up",
x2==1 & x3==0 & x9==1 & x10==-1 ~ "Quintic Down",
x2==1 & x3==0 & x9==-1 & x10==1 ~ "Quintic Down",
x2==0 & x3==1 & x9==-1 & x10==-1 & x8>0 ~ "Quadratic Up Right",
x2==0 & x3==1 & x9==1 & x10==1 & x8>0~ "Quadratic Up Right",
x2==0 & x3==1 & x9==1 & x10==-1 & x8<0 ~ "Quadratic Up Left",
x2==0 & x3==1 & x9==-1 & x10==1 & x8<0~ "Quadratic Up Left",
x2==0 & x3==1 & x9==-1 & x10==1 & x8>0~ "Quadratic Down Right",
x2==0 & x3==1 & x9==1 & x10==-1 & x8>0~ "Quadratic Down Right",
x2==0 & x3==1 & x9==-1 & x10==-1 & x8<0~ "Quadratic Down Left",
x2==0 & x3==1 & x9==1 & x10==1 & x8<0~ "Quadratic Down Left",
)) %>%
mutate( x11 = ifelse(mycolor=="Quadratic Down Right", 1, 0),
x12 = ifelse(mycolor=="Quadratic Up Right", 1, 0),
x13 = ifelse(mycolor=="Cubic Down", 1, 0),
x14 = ifelse(mycolor=="Cubic Up", 1, 0),
x15 = ifelse(mycolor=="Quintic Down", 1, 0),
x16 = ifelse(mycolor=="Quintic Up", 1, 0),
x17 = ifelse(mycolor=="Quadratic Down Left", 1, 0),
x18 = ifelse(mycolor=="Quadratic Up Left", 1, 0))
final.lm <- lm(y ~ x11 + x11: x8 + x11:I( x8^2) + #Quadratic Down Right
x12 + x12: x8 + x12:I( x8^2) + #Quadratic Up Right
x17 + x17: x8 + x17:I( x8^2) + #Quadratic Down Left
x18 + x18: x8 + x18:I( x8^2) + #Quadratic Up Left
x13: x8 + x13:I( x8^3) + #Cubic Down
x14: x8 + x14:I( x8^3) + #Cubic Up
x15: x8 + x15:I( x8^3) + x15:I( x8^5) + #Quintic Down
x16: x8 + x16:I( x8^3) + x16:I( x8^5) #Quintic Up
, data=rdat2)
summary(final.lm)
Call:
lm(formula = y ~ x11 + x11:x8 + x11:I(x8^2) + x12 + x12:x8 +
x12:I(x8^2) + x17 + x17:x8 + x17:I(x8^2) + x18 + x18:x8 +
x18:I(x8^2) + x13:x8 + x13:I(x8^3) + x14:x8 + x14:I(x8^3) +
x15:x8 + x15:I(x8^3) + x15:I(x8^5) + x16:x8 + x16:I(x8^3) +
x16:I(x8^5), data = rdat2)
Residuals:
Min 1Q Median 3Q Max
-12.0607 -2.6280 -0.0454 2.6323 13.0411
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.05800 0.08474 -0.684 0.494
x11 -36.12331 1.72245 -20.972 < 2e-16 ***
x12 35.41434 1.89730 18.666 < 2e-16 ***
x17 -34.54107 2.67490 -12.913 < 2e-16 ***
x18 37.99485 2.35524 16.132 < 2e-16 ***
x11:x8 56.23008 3.19009 17.626 < 2e-16 ***
x11:I(x8^2) -22.77164 1.26613 -17.985 < 2e-16 ***
x8:x12 -52.77655 3.24250 -16.277 < 2e-16 ***
I(x8^2):x12 21.06233 1.21148 17.386 < 2e-16 ***
x8:x17 -49.97923 4.33837 -11.520 < 2e-16 ***
I(x8^2):x17 -19.87236 1.61987 -12.268 < 2e-16 ***
x8:x18 57.71141 3.81713 15.119 < 2e-16 ***
I(x8^2):x18 23.23116 1.38806 16.736 < 2e-16 ***
x8:x13 10.34734 0.52746 19.617 < 2e-16 ***
x13:I(x8^3) -15.12264 0.43847 -34.490 < 2e-16 ***
x8:x14 -9.60837 0.51705 -18.583 < 2e-16 ***
I(x8^3):x14 14.68082 0.41038 35.774 < 2e-16 ***
x8:x15 -2.64218 0.50620 -5.220 1.94e-07 ***
I(x8^3):x15 2.85507 0.30207 9.452 < 2e-16 ***
x15:I(x8^5) -0.45953 0.03993 -11.509 < 2e-16 ***
x8:x16 3.19876 0.50060 6.390 1.96e-10 ***
I(x8^3):x16 -2.92284 0.28742 -10.169 < 2e-16 ***
I(x8^5):x16 0.45776 0.03626 12.625 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.934 on 2577 degrees of freedom
Multiple R-squared: 0.7723, Adjusted R-squared: 0.7704
F-statistic: 397.3 on 22 and 2577 DF, p-value: < 2.2e-16Math Model
\[ Y_i = \beta_{0} + \beta_{1} x_{11} + \beta_{2} x_{12} + \beta_{3} x_{17} + \beta_{4} x_{18} + \beta_{5} x_{11} x_{8} + \beta_{6} x_{11} x_{8^2} + \beta_{7} x_{8} x_{12} + \beta_{8} x_{8^2} x_{12} + \beta_{9} x_{8} x_{17} + \beta_{10} x_{8^2} x_{17} + \beta_{11} x_{8} x_{18} + \beta_{12} x_{8^2} x_{18} + \beta_{13} x_{8} x_{13} + \beta_{14} x_{13} x_{8^3} + \beta_{15} x_{8} x_{14} + \beta_{16} x_{8^3} x_{14} + \beta_{17} x_{8} x_{15} + \beta_{18} x_{8^3} x_{15} + \beta_{19} x_{15} x_{8^5} + \beta_{20} x_{8} x_{16} + \beta_{21} x_{8^3} x_{16} + \beta_{22} x_{8^5} x_{16} + \epsilon_i \]
Plot
palette(c("skyblue","orange","green","purple","steelblue","red","green3","black"))
plot(y ~ x8, data=rdat2, col=as.factor(mycolor))
points(final.lm$fit ~ x8, data=rdat2, col=as.factor(mycolor), pch=16, cex=0.5)
b <- coef(final.lm)
drawit <- function( x11=0, x12=0, x13=0, x14=0, x15=0, x16=0, x17=0, x18=0, i=1){
curve(b[1] + b[2]* x11 + b[3]* x12 + b[4]* x17 + b[5]* x18 + b[6]* x11* x8 + b[7]* x11* x8^2 + b[8]* x8* x12 + b[9]* x8^2* x12 + b[10]* x8* x17 + b[11]* x8^2* x17 + b[12]* x8* x18 + b[13]* x8^2* x18 + b[14]* x8* x13 + b[15]* x13* x8^3 + b[16]* x8* x14 + b[17]* x8^3* x14 + b[18]* x8* x15 + b[19]* x8^3* x15 + b[20]* x15* x8^5 + b[21]* x8* x16 + b[22]* x8^3* x16 + b[23]* x8^5* x16, add=TRUE, xname="x8", col=palette()[i])
}
drawit(1,0,0,0,0,0,0,0,4)
drawit(0,1,0,0,0,0,0,0,6)
drawit(0,0,1,0,0,0,0,0,1)
drawit(0,0,0,1,0,0,0,0,2)
drawit(0,0,0,0,1,0,0,0,7)
drawit(0,0,0,0,0,1,0,0,8)
drawit(0,0,0,0,0,0,1,0,3)
drawit(0,0,0,0,0,0,0,1,5)
with(rdat2, levels(interaction( x11, x12, x13, x14, x15, x16, x17, x18, drop=TRUE)))## [1] "1.0.0.0.0.0.0.0" "0.1.0.0.0.0.0.0" "0.0.1.0.0.0.0.0" "0.0.0.1.0.0.0.0"
## [5] "0.0.0.0.1.0.0.0" "0.0.0.0.0.1.0.0" "0.0.0.0.0.0.1.0" "0.0.0.0.0.0.0.1"